subroutine backup_hydro(filename, filename_desc)
  use amr_commons
  use hydro_commons
  use dump_utils, only : dump_header_info, generic_dump, dim_keys
  use mpi_mod
  implicit none
#ifndef WITHOUTMPI
  integer :: dummy_io, info2
#endif

  character(len=80), intent(in) :: filename, filename_desc

  integer :: i, ivar, ncache, ind, ilevel, igrid, iskip, istart, ibound
  integer :: unit_out, unit_info
  real(dp) :: d
#ifdef SOLVERmhd
  real(dp) :: A, B, C
#endif
  integer, allocatable, dimension(:) :: ind_grid
  real(dp), allocatable, dimension(:) :: xdp
  character(LEN = 5) :: nchar
  character(LEN = 80) :: fileloc
  integer, parameter :: tag = 1121
#if NENER > 0
  integer :: irad
#endif
  logical :: dump_info_flag
  integer :: info_var_count
  character(len=100) :: field_name

  if (verbose) write(*,*)'Entering backup_hydro'

  call title(myid, nchar)
  fileloc = TRIM(filename)//TRIM(nchar)

  ! Wait for the token
#ifndef WITHOUTMPI
  if (IOGROUPSIZE > 0) then
     if (mod(myid-1, IOGROUPSIZE) /= 0) then
        call MPI_RECV(dummy_io, 1, MPI_INTEGER, myid-1-1, tag,&
             & MPI_COMM_WORLD, MPI_STATUS_IGNORE, info2)
     end if
  end if
#endif

  open(newunit=unit_out, file=fileloc, form='unformatted')

  if (myid == 1) then
     open(newunit=unit_info, file=filename_desc, form='formatted')
     call dump_header_info(unit_info)
     info_var_count = 1
     dump_info_flag = .true.
  else
     dump_info_flag = .false.
  end if

  write(unit_out) ncpu
  if(strict_equilibrium>0)then
     write(unit_out) nvar_all+2
  else
     write(unit_out) nvar_all
  endif
  write(unit_out) ndim
  write(unit_out) nlevelmax
  write(unit_out) nboundary
  write(unit_out) gamma
  do ilevel = 1, nlevelmax
     do ibound = 1, nboundary+ncpu
        if (ibound <= ncpu) then
           ncache = numbl(ibound, ilevel)
           istart = headl(ibound, ilevel)
        else
           ncache = numbb(ibound-ncpu, ilevel)
           istart = headb(ibound-ncpu, ilevel)
        end if
        write(unit_out) ilevel
        write(unit_out) ncache
        if (ncache > 0) then
           allocate(ind_grid(1:ncache), xdp(1:ncache))
           ! Loop over level grids
           igrid = istart
           do i = 1, ncache
              ind_grid(i) = igrid
              igrid = next(igrid)
           end do
           ! Loop over cells
           do ind = 1, twotondim
              iskip = ncoarse+(ind-1)*ngridmax
              do ivar = 1, neul-1
                 if (ivar == 1) then
                    ! Write density
                    do i = 1, ncache
                       xdp(i) = uold(ind_grid(i)+iskip, 1)
                    end do
                    field_name = 'density'
                 else
                    ! Write velocity field
                    do i = 1, ncache
                       xdp(i) = uold(ind_grid(i)+iskip, ivar)/max(uold(ind_grid(i)+iskip, 1), smallr)
                    end do
                    field_name = 'velocity_' // dim_keys(ivar - 1)
                 end if
                 call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
              end do
#ifdef SOLVERmhd
              do ivar = 6, 8 ! Write left B field
                 do i = 1, ncache
                    xdp(i) = uold(ind_grid(i)+iskip, ivar)
                 end do
                 field_name = 'B_' // dim_keys(ivar - 6 + 1) // '_left'
                 call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
              end do
              do ivar = nvar+1, nvar+3 ! Write right B field
                 do i = 1, ncache
                    xdp(i) = uold(ind_grid(i)+iskip, ivar)
                 end do
                 field_name = 'B_' // dim_keys(ivar - (nvar+1) + 1) // '_right'
                 call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
              end do
#endif
#if NENER > 0
              ! Write non-thermal pressures
              do ivar = nhydro+1, nhydro+nener
                 do i = 1, ncache
                    xdp(i) = (gamma_rad(ivar-nhydro)-1d0)*uold(ind_grid(i)+iskip, ivar)
                 end do
                 write(field_name, '("non_thermal_pressure_", i0.2)') ivar-nhydro
                 call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
              end do
#endif
              ! Write thermal pressure
              do i = 1, ncache
                 d = max(uold(ind_grid(i)+iskip, 1), smallr)
                 xdp(i) = uold(ind_grid(i)+iskip, neul)
                 xdp(i) = xdp(i)-0.5d0*uold(ind_grid(i)+iskip, 2)**2/d
#if NDIM > 1 || SOLVERmhd
                 xdp(i) = xdp(i)-0.5d0*uold(ind_grid(i)+iskip, 3)**2/d
#endif
#if NDIM > 2 || SOLVERmhd
                 xdp(i) = xdp(i)-0.5d0*uold(ind_grid(i)+iskip, 4)**2/d
#endif
#ifdef SOLVERmhd
                 A = 0.5d0*(uold(ind_grid(i)+iskip, 6)+uold(ind_grid(i)+iskip, nvar+1))
                 B = 0.5d0*(uold(ind_grid(i)+iskip, 7)+uold(ind_grid(i)+iskip, nvar+2))
                 C = 0.5d0*(uold(ind_grid(i)+iskip, 8)+uold(ind_grid(i)+iskip, nvar+3))
                 xdp(i) = xdp(i) - 0.5*(A**2+B**2+C**2)
#endif
#if NENER > 0
                 do irad = 1, nener
                    xdp(i) = xdp(i)-uold(ind_grid(i)+iskip, nhydro+irad)
                 end do
#endif
                 xdp(i) = (gamma-1d0)*xdp(i)
              end do
              field_name = 'pressure'
              call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
#if NVAR > NHYDRO+NENER
              ! Write passive scalars if any
              do ivar = nhydro+1+nener, nvar
                 do i = 1, ncache
                    xdp(i) = uold(ind_grid(i)+iskip, ivar)/max(uold(ind_grid(i)+iskip, 1), smallr)
                 end do
                 if (metal .and. imetal == ivar) then
                    field_name = 'metallicity'
                 else
                    write(field_name, '("scalar_", i0.2)') ivar - nhydro - 1 - nener
                 end if
                 call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
              end do
#endif
              if(strict_equilibrium>0)then
                 do i = 1, ncache
                    xdp(i) = rho_eq(ind_grid(i)+iskip)
                 end do
                 field_name = 'equilibrium_density'
                 call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
                 do i = 1, ncache
                    xdp(i) = p_eq(ind_grid(i)+iskip)
                 end do
                 field_name = 'equilibrium_pressure'
                 call generic_dump(field_name, info_var_count, xdp, unit_out, dump_info_flag, unit_info)
              endif
              ! We did one output, deactivate dumping of variables
              dump_info_flag = .false.
           end do
           deallocate(ind_grid, xdp)

        end if
     end do
  end do
  close(unit_out)

  if (myid == 1) close(unit_info)
  ! Send the token
#ifndef WITHOUTMPI
  if (IOGROUPSIZE > 0) then
     if (mod(myid, IOGROUPSIZE) /= 0 .and.(myid .lt. ncpu)) then
        dummy_io = 1
        call MPI_SEND(dummy_io, 1, MPI_INTEGER, myid-1+1, tag, &
             & MPI_COMM_WORLD, info2)
     end if
  end if
#endif


end subroutine backup_hydro
